Nonlinear terahertz emission in semiconductor microcavities 
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We consider the nonlinear terahertz emission by the system of cavity polaritons in the regime 
of polariton lasing. To account for the quantum nature of terahertz- polariton coupling we use 
Lindblad master equation approach and demonstrate that quantum microcavities reveal rich variety 
of the nonlinear phenomena in terahertz range, including bistability, short THz pulse generation 
and THz switching. 
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Introduction. THz band remains the last region of elec- 
tromagnetic spectrum which does not have wide appli- 
cation in modern technology due to lack of solid state 
source of THz radiation which is compact, reliable and 
scalable Fundamental objection preventing realiza- 
tion of such source is small rate of spontaneous emission 
of the THz photons. According to Fermi Golden rule this 
rate is about tens of inverse milliseconds, while lifetime of 
the photoexcited carrier typically lies in picosecond range 
due to the efficient interaction with phonons 0, H[ • Spon- 
taneous emission rate can be increased by application of 
Purcell effect when emitter of THz is placed in cavity for 
THz mode 0, but even in this case cryogenic tem- 
peratures are required to provide quantum efficiency of 
the order about one percent for typical quantum cascade 
structure. 

Recently it was proposed that the rate of spontaneous 
emission for THz photons can be additionally increased 
by bosonic stimulation if radiative transition occurs into 
a condensate state Q. One example is a transition be- 
tween upper and lower polariton branches in semiconduc- 
tor microcavity in the regime polariton lasing. Unfortu- 
nately, the radiative transition accompanied by emission 
of THz photon between upper and lower polariton modes 
is forbidden, since these states have the same parity. 
Nevertheless, such transition becomes possible if upper 
polariton state is mixed with exciton state of different 
parity. Amplification of spontaneous emission by Pur- 
cell effect together with bosonic stimulation increase the 
rate of spontaneous emission by several orders of magni- 
tude, making it comparable with the rate of scattering 
with acoustic phonon. Consequently, effective emission 
of THz radiation can occur [f| . 

It is well known that strong polariton-polariton inter- 
actions in microcavities make it possible to observe pro- 
nounced nonlinear effects for the intensities of the pump 
orders of magnitude smaller than in other nonlinear opti- 
cal systems. Among them arepolariton superfluidity [J , 
bistability and multistability ja, 0] , soliton formation [lC| 
and others. One can expect that polariton-polariton in- 



teractions will as well strongly affect the process of THz 
emission. The quasiclassical approach based on Boltz- 
mann equations, used in Ref@ cannot provide a correct 
account of a coherent interaction of THz photons and po- 
laritons, and can not be used for satisfactory description 
of nonlincarities in the considered system. The develop- 
ment of more exact quantum formalism is thus needed. 
This paper is aimed at building such a formalism, which 
accounts for the following physical processes: coherent 
polariton-THz photon interaction, polariton- polariton 
interaction leading to the blucshift of the polariton modes 
and coupling of the polaritons with acoustic phonons. 
The development of such description is timely in light 
of intensive studies of ultrastrong light-matter coupling 
Tllfl2 |. single cycle THz generation [3, EH, mtersub- 



band cavity polariton physics HH and control of the 
phase of THz radiation in both inorganic and organic 
structures |l£ 
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Formalism. We consider a model system consisting of 
a lower polariton state with the energy 6l, upper hybrid 
state with the energy eu, THz cavity mode with the en- 
ergy ct and incoherent polariton reservoir coupled with 
upper and lower polariton states via phonon-assisted pro- 
cess (see Fig[T|). 

The Hamiltonial of the system written in terms of the 
operators of secondary quantization for upper polaritons 
(au,a>u)> l° wer polaritons (a£,a£), TH photons (c, c + ), 
reservoir states (a,Rk, a^/) and acoustic phonons bk, b^ 
can be represented as a sum of four terms: 



H = H + H, 



pol—pol 



Ht + He 



(i) 



The first term 

H = t L a\a L + eua+ajj + e T c + c + ^ £ka^ k a flk (2) 



corresponds to the energy of uncoupled upper and lower 
polaritonic states, THz mode and polariton reservoir. 
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where and 6k denote operators of creation and anni- 
hilation of phonons with wavevector k, D\ and D2 are 
the polariton-phonon interaction constants. 

Keeping in mind that interactions described by 
Hq, Hp i_ po i, Ht are of coherent nature, while phonon 
assisted interactions (Hr) with reservoir destroy coher- 
ences the dynamics of the density matrix of system p 
is described by Lindblad master equation, analogical to 
those obtained in Refs 2(| HI (see also supplementary ma- 
terial). 



^ = ^[p;H + H V oi- V oi + H T ] + (6) 



FIG. 1: The scheme of transitions in THz emitting cavity. 
The upper polariton is mixed with dark exciton state due to 
the application of the gate voltage V a . The radiative transi- 
tion between the upper hybrid state \U) and lower polariton 
state \L) thus becomes possible. Upper hybrid and lower po- 
lariton states are also coupled with an incoherent reservoir of 
the polaritons via phonon-assisted process. 



The second term 

Hpoi-poi = U LLO.\a\ahaL + Uuuo-u a u a u a u +(3) 

+2U UL a u a u a L ~a L + 2J {U L Ra[a L + U UR a^}au) a^ k a# k 
k 

describes polariton- polariton interaction. The intcrac- 



5 -^L{2{H+pH- R +H-pHt) 



{H + R H R + H R H+) p - p (H+Hr + H R H+) } 
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tion constants can be estimated as U;, 



X?X?U, where 



Xi are Hophcld coefficients giving the percentage of the 
exciton fraction in the polariton states. Xjj and Xl are 
determined by cavity geometry, and we took Xr = 1 sup- 
posing that the reservoir is purely excitonic. The matrix 
element of the exciton- exciton scattering can be esti- 
mated as U re 6Esa B /S with Eb and cle being the exci- 
ton binding energy and Bohr radius respectively, and S 
the area of the system . 
The third term 



where La is Linblad operator defined by the formula 
La = 2ApA~ — A + Ap — pA + A and tl, tu, tr and r 
are lifetimes of lower polaritons, upper polaritons, po- 
laritons in the reservoir and THz photons, and P and I 
are pumping intensities of upper polariton state and ter- 
ahertz mode. The delta function 8&e denotes the conser- 
vation of energy in the process of phonon scattering. The 
first line accounts for the coherent processes in the sys- 
tem, the second and third lines correspond to the phonon- 
assisted coupling with incoherent reservoir of the polari- 
tons, the last line accounts for the pump and the decay. 

The equations for the populations of polariton states 
and terahertz photons can be obtained as 



Tr 



dp 
'~dt 



(7) 



Hn 



V T (ata L c + a u a L 'c + ) 



(4) 



describes radiative THz transion between upper and 
lower polariton states. The matrix element of the THz 
emission can be estimated using a standard formula for 
the coupling constant of the dipole transition with con- 
fined electromagnetic mode, Vt = uP'd^J hn/2ir 3 eoC 3 , 
where d is matrix element of the radiative transition 
and n- refrective index of the terahertz cavity (See e.g. 
ReiJH). 

Interaction between upper and lower polariton states 
and incoherent reservoir is described by the fourth term: 

Hr = H+ + Hx = Di ]T {a v a% k b+ + a+ajoA) + (5) 



Using the mean field approximation, one gets the 
closed system of the dynamic equations for the occu- 
pancies n L = (a£a L ), n v = (a^au),n R u = (a^ k a flk ) 
and n = (c + c) connected by the correlators oclu = 
(a~^auc + ),auL = (aifljc) = a* LU (sec supplementary 
material) 

d t n L = -2^Im (a UL ) - ^ + (8) 

n tl 
+W 2 ^{(n £ + 1) "flk (ni h + l)-n L (n Rk + 1) nf }; 



cW=2^Im(o m )-^ + P+ (9) 



aLaR k D k 



-W x Y,i( n u + 1) «flk< /l - n v (n Rk + 1) (nf + l) }; 
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dtnRq = ———^+ (10) between upper and lower polariton modes equal to 16 

Tr meV (which corresponds to 4 THz) and embedded into 

Wi ^2{nu (n Rk + 1) (n^ h + l \ - (riu + 1) n Rk n^} + THz cavity with eigen frequency slightly different from 

k Qr and having a quality factor Q = 100 U 23 1. Let 



, Jlr \ ■> r i , 1 % P h / ,i / »ii , ,\, us assume that initially the system is characterized by 

+W 2 > {n L (n Ris . + 1) n£ - (n L + l)n Rk n£ + 1 }; p i n i to w ,,,, 

V / zero population of polaritons and 1Hz photons. When 

the constant non-resonant pump of the upper polariton 

state is switched on, the number of THz photons n starts 

V T n to increase until it reaches some equilibrium level defined 

= ~ 2_ fi" Im ^ auL ' ~ ~ + 1 ^ by the radiative decay of polaritons and escape of THz 

radiation from the cavity, as it is shown at FigPJ] 



d t a UL = - (eu - e L - e T ) a UL h (12) 

v T 

+i-r-{{nu + 1) n L n - n v (n L + 1) (n + 1)} + 
+{Wi (-nia* - K h - l) + W 2 £ (n Rk - nl h ) }a UL . 

k k 

In the above expressions r~ l r = 1 + 1 + t^ 1 + t _1 , 
Vt ~ 1 fJ-eV is a coupling constant between polaritons 
and terahertz photons and W± t 2 ~ 2ps _1 are transition 
rates between the reservoir and upper/lower polariton 
states determined by polariton-phonon interaction con- 
stants, Wi t 2 ~ |-Di,2| 2 - Note, that chatracteristic time of 
tcrzhetz photon emission is about three orders of mag- 
nitude smaller than characteristic time of the scattering 
with acoustic phonons. However, THz emission is dra- 
matically inhanced by bosonic stimulation and becomes 
dominant mechanism for sufficiently strong pumps, n^ 1 
gives the occupancies of the phonon mode determined 
by Bose distribution function. For simplicity of the cal- 
culations in the present paper we consider the reservoir 
to consist from N identical states (N = 3 • 10 5 ). Note, 
that if coherent interaction is switched off by equating 
dctjjhldt = the system of the equations we use trans- 
forms into the system of Boltzmann equations considered 
in RefE 

The renormalized energies of the upper and lower po- 
lariton states are determined by their blueshifts arising 
from polariton-polariton interactions and read 

/ \ 

% = e v + 2 Uuunu + U UL n L + U UR ^ n R^ > ( 13 ) 

e L = e L + 2 I U LL n L + U UL nu + U LR ^ ™Rk J • (14) 

Due to the difference of the Hopfield coefficients for the 
upper and lower polariton states, the difference eu — ii 
depends on the polariton concentrations and thus is de- 
termined by the intensity of the pump P. This depen- 
dence can have important consequencies, allowing for the 
onset of the bistability in the system (see below). 

Results and discussions. We consider a planar GaAs 
microcavity in strong coupling regime with Rabi splitting 
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FIG. 2: Time evolution of terahertz photons number at zero 
temperature for different pumps: 4500 ps _1 (red), 5000 ps _1 
(green), 5500 ps _1 (blue) and 6000 ps _1 (black). Inset shows 
evolution of THz photons number for the constant pump P = 
6 10 3 ps _1 for different temperatures: 1 K (red), 10 K (green) 
and 20 K (blue). 

Equilibrium value of the THz population n as a func- 
tion of pumping P demonstrates threshold-like behavior. 
For high enough temperatures, below the threshold the 
dependence of n on P is very weak. When pumping 
reaches the certain threshold value, polariton conden- 
sate is formed in the lower polariton state, radiative THz 
transition is amplified by bosonic stimulation, and the 
occupancy of THz mode increases superlinearly together 
with the occupancy of lower polariton state (Fig^l 
blue curve). This behavior is qualitatively the same as 
in the approach operating with semiclassical Boltzmann 
equations. However, the decrease of temperature leads 
to the onset of the bistability and hysteresis in the de- 
pendence n(P). The bistable jump occurs when the in- 
tensity of the pump tunes eu — £l into the resonance 
with the cavity mode ex- The parameters of the hystere- 
sis loop strongly depend on the temperature (FigJ3|). It 
is very pronounced and broad for low temperatures, nar- 
rows with the increase of the temperature and disappears 
completely at T sa 20K. 

Coherent nature of the interaction beween excitons and 
THz photons makes possible the periodic exchange of 
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FIG. 3: Dependence of occupancy of the THz mode on pump 
in equilibrium state for different temperatures: 1 K (red), 10 
K (green) and 20 K (blue). 

the energy between polaritonic and photonic modes and 
oscillatory dependence of the THz signal in time. Fig|3] 
shows temporal evolution of the occupancy of THz mode 
after excitation of the upper polariton state by a short 
pulse having a duration of about 2 ps. It is seen that the 
occupancy of THz mode reveals a sequence of the short 
pulses having duration of dozens of ps with amplitude 
decaying in time due to escape of THz photons from a 
cavity and radiative decay of polaritons. The period of 
the oscillations is sensitive to the number of the injected 
polaritons No and decreases with increasing of No. If 
the lifetime of polaritons is less than the period of the 
oscillations, single pulse behaviour can be observed as it 
is shown in the inset of FigJH Appropriate choice of the 
parameters can ultimately lead to a generation of THz 
wavelets composed of one or several THz cycles, which 
makes polariton-THz system suitable for application in 
a sort pulse THZ spectroscopy. 

If the system of coupled THz photons and cavity po- 
laritons is in the state corresponding to the lower branch 
of the S-shaped curve in the bistability region, illumina- 
tion of the system by a short THz impulse I{t) can induce 
its switching to the upper branch, as it is demonstrated 
at Figf5] One sees, that the response of the system is 
qualitatively different for different values of the pump P. 
If P lies outside the bistablity region, the application of 
THz pulse leads to a short increase of the n but subse- 
quently the system relaxes to its original state (red and 
blue curves). However, when the system is in the bista- 
bility regime the switching occurs. Note that this effect 
is of quantum nature and cannot be describes using the 
approach based on semiclassical Boltzmann equations de- 
veloped in RefE 

Conclusion. We considered the system of coupled cav- 
ity polaritons and THz photons using the approach based 
on generalized Lindblad equation for the density ma- 
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after 

the arrival of the short excitation pulse (dashed black line). 

FIG. 4: The temporal dependence of the terahertz mode oc- 
cupancy. The background pump is switched off, lifetimes 
tl = tu = 50 ps. The temperatures are: 1 K (red line), 
10 K (green) and 20 K (blue). Inset: T=1K, different life- 
times of the polariton states: 15 ps (red), 20 ps(green) and 
25 ps(blue). 
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FIG. 5: System response on a short single impuls (FWHM= 
40 ps, black/dotted curve) for different values of the back- 
ground pump: 4600 ps -1 (red), 5750 ps -1 (green), 6500 ps -1 
(blue). Switching occurs only when the background pump 
corresponds to the bistable region (compare with Fig. 2). 

trix. We showed that such system demonstrates a va- 
riety of intriguing nonlinear effects, including bistability, 
THz swithching and generation of short THz wavelets. 
The work was supported by Rannis "Center of excel- 
lence in polaritonics" , FP7 IRSES "POLAPHEN" and 
"POLALAS" projects, Russian Fund of Basic Research 
and COST "POLATOM" program. 
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SUPPLEMENTARY MATERIALS 



In the present supplementary appendix we present a derivation of the quantum kinetic equations for the system of 
cavity polaritons coupled with a terahertz (THz) cavity mode based on the Lindblad approach for the density matrix 
dynamics. The method we develop is general and can in principle be applied to any system of interacting bosons in 
contact with a phonon reservoir, for example, a polaritonic channel [18,19] or a condensate of indirect excitons. 



The Lindblad approach 



The system of polaritons, phonons and THz cavity photons is described by its density matrix p, for which we apply 
Born approximation factorizing it into the phonon part which is supposed to be time-independent and corresponds to 

the thermal distribution of acoustic phonons p p h = exp |— j3H p h j, and the part describing polaritons and THz cavity 
photons pcav whose time dependence should be determined, p = p p h <8> p ca v ■ Our aim is to find dynamic equations 
for the time evolution of the occupancies of the upper and lower polariton states and the THz cavity mode: 

n i {t) = Tr\a\a i p{t)) = {a\a t ) (15) 

where aj and cij are operators of the upper and lower polaritons (i = U and i = L respectively) and THz cavity photons 
(i = T). In our consideration we neglect spin of the cavity polaritons since out goal is to find the effects of bistability 
and switching in the THz emitter and spin degree of freedom is not expected to introduce any qualitatively new physics 
from this point of view. It should be noted, however, that introduction of spin into the model is straightforward. 
The total Hamiltonian of the system can be represented as a sum of two parts 

H = hW+H& (16) 

where the first term HW describes the "coherent" part of the evolution, corresponding to free polaritons, cavity 
photons and polariton-polariton interactions, and the second term corresponds to the dissipative interaction 

with acoustic phonons. The two terms affect the polariton density matrix in a qualitatively different way. The effect 
of the coherent part on the evolution of the density matrix is described by the Liouville-von Neumann equation 

i{d tP ) (1) = [ffW;?] (17) 

Polariton-phonon scattering corresponds to the interaction of the quantum polariton system with a classical phonon 
reservoir. It is of dissipative nature, and thus straightforward application of the Liouville-von Neumann equation is 
impossible. Introduction of dissipation into quantum systems is an old problem, for which there is still no single well 
established solution. In the domain of quantum optics, however, there are standard methods based on the Master 
Equation techniques. In the following we give a brief outline of this approach applied for a dissipative polariton 
system. 
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The Hamiltonian of interaction of polaritons with acoustic phonons in Dirac representation can be represented as 

H^(t) = H-(t)+H+{t) = DuY,e i{Eu - Enh)t (aua+ k +a+a Rk ) {b k e~ 1 ^ 1 + 6+e^ 4 ) + (18) 

fc 

+D L Y,j {ERk - EL)t (a L a+ k +ata Rk ) foe"*"** + b+e^ 



where a, are the operators for polaritons, b k are the operators for phonons, Ei and w k are the dispersion relations of 
polaritons and acoustic phonons respectively, Di are the polariton-phonon coupling constants. In the last equality we 
separated the terms H + where a phonon is created, containing the operators 6 + , from the terms H~ in which it is 
destroyed, containing operators b. 

Now, one can consider a hypothetical situation when polariton-polariton interactions are absent, and all redistri- 
butions of the polaritons are due to the scattering with a thermal reservoir of acoustic phonons. One can rewrite the 
Liouville-von Neumann equation in an intcgro-diffcrcntial form and apply the so-called Markovian approximation 



(dtp) 



(2) 



1 



(19) 



Sae 2{H+pH- +H-pH~ 



H + H~ +H-H+) p-p(H + R- + H~H 



Lh^Pi 



where in the last line symbol L H m denotes a Lindblad dissipative operator corresponding to the Hamiltonian H^ 2 \ 
The coefficient Sae corresponds to the energy conservation and has dimensionality of inverse energy divided by square 
of a Plank constant. In calculations we estimate Sae as being proportional to the inverse broadening of the polaritons 
states (as it is usually done in calculation of the transition rates in semiclassical Boltzmann equations using Fermi 
golden rule). For time evolution of the mean value of any arbitrary operator (A) = Tr(pA) due to scattering with 
phonons one thus has: 



{d t (A)} {2) = Sae (<[#"; [A; £+]]) + <[#+; $£-]])) 



(20) 



Putting A = afcii in this equation we get the contributions to the dynamic equations for the occupancies coming 
from polariton-phonon interactions, which are nothing more than the standard semi-classical Boltzmann equations 
describing the thermalization of a polariton system. 

Combined together, coherent and incoherent contributions result in a following master equation for the density 
matrix: 



(21) 



Coherent part 

Let us consider the coherent part. Hamiltonian here = Hq + Ht + H po i^ po i (let us omit symbols " ~" over 
operators hereafter.), where 



is the free-polariton Hamiltonian, 



Ho = eia^aL + eira^au + £tc + c 



(22) 



H T = V T (a^a L c + aua^c + ) 



(23) 



is the polaritons-to-THz photons interaction term and 
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Hp i-p i = V LL a\a\a L a L + Uuua^a^auau + 2U UL au a ua\ j a L + 2U UR aj J a u a\ lk a Rk + 2U LR a^a L a^ k a Rk «(24) 

k k 

w 2 [U LL {a\ j a L {ala L )) + U vv (a^auia^au)) + U UL (a^aLia^au)) + U UL (a^au (a£a L )) + 
+U UR (a^}a u (a R ' k a R k)) + U UR (a^ k a Rk (a^au)) + U LR (a\ j a L {a\ fk a Rk )) + U LR (a+ k a Rk (ala L ))] « 

w 2 (u LL n L + U UL nu + ^ U LR n Rk ^j a\a h + 2 (ujjijnu + U UL n L + ^ U UR n Rk J a^a v + 
+2 (U UR nu + U LR n L ) ^ a R k a Rk = U\a\a L + U 2 a^a u + U 3 2J a Rk a Rk ; 

k k 

is the polariton-polariton scattering Hamiltonian. Coming from the first to the second lines of this expression we 
used mean-field approximation. Here we also introduced three new coefficients: 



Ux = 2 ^Ulltil + Uultiu + ^2 ULRn R kj , 

/] UuRn R k , 



1) Lower and upper polaritons occupancies riL,nu, 
Consider the dynamic equation for til as an example. One has 

h(d t n L ) {1) =iTY{a+a L [p-H^]}- 

h{d t n L ) Q = iTr{a+a L [p; H }} = iTr{p[e L a+a L +e u a^a u + e T c + c]a[a L ]} = 0; 
h(d t n L ) T = iV T TT{p[a^a L c + a v a^c + ;a\a L ]} = iV T (a UL - a^} L ) = -2V T lm(a UL ) ; 

ft (dtn L ) pol _ pol = iTT{p[Uia+a L + U 2 a^au + U 3 ^ a^ k a Rk ; aJa L ]} = 0. 



Finally, 



(25) 



( Uuunu + U UL n L + ^ U UR n Rk j , (26) 
U 3 = 2(U UR nu + U LR n L ). (27) 



ft (d t n L ) {1) = -2V T lm {a UL ) . (28) 



The equation for the upper polariton occupancy is obtained in a similar way. 

2) Reservoir n Rk , 

h(d t n Rk ) (1) = iTT{a+ k a Rk [p;H^}} = 0. (29) 

3) Terahertz cavity occupancy rij, 

h{d f n T ) {1) =tTr{c+c[p-H^}}; 

ft(d t n T ) o = ;Tr{ C + c[p;ir o ]} = 0; 
h{d t n T ) T = iVT^{p[a^aijc + a v a^;c + c\] = iV T (a UL - a^j L ) = ~2V T lm(a UL ) ; 

ft {dtn T ) pol _ pol = iTr{p[C/ia+a L + U 2 a^au + U 3 ^ a^ k a Rk ; c + c}} = 0. 

fc 

Finally, 

ft (d t n T ) (1} = -2V T Im (a UL ) . (30) 
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4) Correlators ajjL, 

h(d t a UL ) {1) = iTr{a+a L c[p;H^]}; 

h(d t a UL ) = iTv{a^a L c[p; H Q ]} = iTr{p[e L a L a L + eua^au + £tc + c; a^a L c]} = i {-e L + e v - e T ) a UL ; 
h(d t a UL ) T = iV T Tr{p[aua L c + aua L c + ;a^a L c]} = iV T {(a u a L a^ r a L c r c) - {a~}}aLauat cc+ )) = 

= -iVr{(nu + 1) n L n T - n v (n L + 1) (n T + 1)}; 
ft {d t a UL ) pol _ pol = iTv{p[Uia L a L + U 2 af J a u + U 3 ^ a+ k a Rk ; a^a L c}} = 

k 

= iUiTr{p[a L a L ] a^a L c]} + iU2^{[a^au; a^}a L c]} = iUi (-a UL ) + iU 2 a UL - 

Finally, 

ft (dtauL)^ 1) = i( e u - £l - £t) a UL + 2V T {{nu + 1) n L n T - (31) 
nu (n L + 1) (n T + 1)} + i{2U uu n u + U UL (n L - nu) - 2U LL n L + (U UR - U LR ) ^ n Rk }a UL ; 

k 

Decoherent part 

To get explicit expresions for the dynamics of ul, nu, n Rk , nx and ajjL due to decoherent processes of interaction 
with the reservoir let us consider Liouville-von Neumann equation for the density matrix after the Born-Markov 
approximation Eq. (6) and the simplest case when only three states L, U and R k are present. 

In this case, leaving energy- conserving terms only one gets 

H + = D x J2 a ua+ R b+ + D 2 J2 ata kR b+ (32) 

k k 

H~ = Pi a+a kR b k + D 2 ^ a L a kR h (33) 

k k 

The application of Eq. (6) gives the following results: 

1) Lower and upper branch polariton occupancies n^ny, 
For til one has: 

(d t n L ) {2) = S AE ([#"; [a+a L ;H+]} + [H+; [a+a L ; H']}) = 2S AE [H+; [a+a L ; #"]]; 
[a L a L ; H~] = [a L a L ; Di 2J a+a Rk b k + D 2 ^ a L a^ k bk] = -D 2 ^ aLa+ fe 6 fc ; 

k k k 

[H+; [a+a L ;H-]} = ^ a v a+ k b+ + D 2 ^ a+a Rk b+; -D 2 ^ a L a+ k b k ] = 
\ D 2 a L a Rk b+; -D 2 ^ a L a+ k b k ] = D\ ^{(aJa L + l) a+ k a Rk (b+b k + l) - a+a L (a^ k a Rk + l) b+b k }; 

Finally, 

h 2 (d t n L )W = w 2 ]T{(n L + 1) («f + l) - nz (n Rfc + 1) nf }, (34) 

k 

where we introduced W 2 = 28 AE D 2 . 



The equation for nu is easily obtained in an analogical way. 
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2) Reservoir n Rk = (a R ~ k a Rk }, 

(d t n Rk ) {2) = 26 AE [H+; [a+ k a Rk ; H~]]; 
[ a Rk a Rk;H~] = [a^ k a Rk ;Di ^ a+,a Rk b k + D 2 2J a L a^ k b k ] = -D 1 ^ a+,a Rk b k + D 2 2J a La Rk b k ; 

lc A,' A? 

[# + [ a j?fc a -R^ ^~]] = aua Rk b k + -°2 ^ a+]a Rk b+; -L>i ^ a^a Rk b k + D 2 ^ a L a Rk b k] ~ 

A.* /c 

« £>i ^{a^ac/ (a+ k a Rk + l) (&+6 fe + l) - (a+a^j + l) a+ k a Rk b+b k } + 

k 

+ D l J2Wl a L ( a Rk a Rk + 1) fofc fofe - (aJaL + l) aflfeORfc (frfc&fc + 1)}, 
fe 

where coming between the third and the fourth lines we neglected the off- diaginal elements of the phonon density 
matrix, supposing (b+b k >) = nf5 k , k > 
Finally, 

(d t n Rk )W « Wi (n Kfe + 1) (nf + l) - (nc; + 1) n Rk nf} + (35) 

+W 2 ^{n L (n Rk + 1) nf - (n L + 1) n Rk (nf + l) }. 

k 

3) Terahertz cavity occupancy ut = (c + c), 

{d t n T f ] = S AE ([H-; [c+c; H+}} + [H+; [c+c; H~]}) = 0. (36) 

4) Decoherent part of the correlator auL = (a^axc), 

[a+ra L c; H~] = 0; 

[a+ra L c; H+] = [a^a L c; D 1 ^ a ua Rk b t + D ^ X! a t a Rkb k ] = - ^ Dia R ~ k a L cbt + ^2 ^ a+ra Rk cb+; 

A? A^ A^ Aj 

[a^a L c; = [D 1 ^ a u a Rkbk + D 2 ^ a La Rk b k r, -D x ^ a Rk a Lcbt + D 2 ^ a+ r a Rk cb+] « 

~ £>f X^"L a Hfe a «fe ( a C/ a i c ) fo fe & fc - { a Rk a Rk + l) a^OLC (6^fe fe + l)} + 

k 

+D\ Y{a+, k a Rk a+ra L c (b+b k + l) - (a+, k a Rk + l) a^a L cb+b k } = 
= -D? / J{-«flfcQ-RfcQqai,c - a^a L c (&^6 fc + l)} + L>2 5Z]{ a /ifc a «fe a c/ a i c _ a+,a L cb+b k }. 

k k 

Finally, 

(dta UL ) {2) = -y- X! _ _ + ~ n fc h ) a£/i ' ( 37 ) 

it fc 

After merging the equations for the coherent processes with the equations for the incoherent phonon-scattering 
processes and adding finite lifetimes, background and impulse pumps (see main text of the Letter) we solve this 
self-consistent set of equations and eventually find the evolu tion of the THz photons occupancy. 



